function matrix_A = assemble_matrix_A_2d(node, elem, function_cor, dN, basis_type, trial_type, test_type, N_gauss_int2d)

matrix_A = sparse(size(node,1), size(node,1));

for n = 1:size(elem,1)
    E = node(elem(n,:),:);
    switch basis_type
        case {"P1", "P1b", "P2"}
            V = E(1:3,1:2);
            phi = function_transform_2d([0,0; 1,0; 0,1], dN, V);
            psi = function_transform_2d([0,0; 1,0; 0,1], dN, V);
        case {"Q1", "Q1b", "Q2"}
            V = E(1:4,1:2);
            phi = function_transform_2d([-1,-1; 1,-1; 1,1; -1,1], dN, V);
            psi = function_transform_2d([-1,-1; 1,-1; 1,1; -1,1], dN, V);
    end
    switch trial_type
        case "dx"
            phi = function_symbolic_computing_2d(phi, "dx");
        case "dxdx"
            phi = function_symbolic_computing_2d(psi, "dxdx");
        case "dy"
            phi = function_symbolic_computing_2d(phi, "dy");
        case "dydy"
            phi = function_symbolic_computing_2d(psi, "dydy");
    end
    switch test_type
        case "dx"
            psi = function_symbolic_computing_2d(psi, "dx");
        case "dxdx"
            psi = function_symbolic_computing_2d(psi, "dxdx");
        case "dy"
            psi = function_symbolic_computing_2d(psi, "dy");
        case "dydy"
            psi = function_symbolic_computing_2d(psi, "dydy");
    end
    switch basis_type
        case {"P1", "P1b", "P2"}
            for alpha = 1:size(elem,2)
                for beta = 1:size(elem,2)
                    f = @(xi,eta) function_cor(xi,eta).*phi{alpha}(xi,eta).*psi{beta}(xi,eta);
                    r = gauss_int2d_tri(f, V(1,:), V(2,:), V(3,:), N_gauss_int2d);
                    i = elem(n,beta);
                    j = elem(n,alpha);
                    matrix_A(i,j) = matrix_A(i,j) + r;
                end
            end
        case {"Q1", "Q1b", "Q2"}
            for alpha = 1:size(elem,2)
                for beta = 1:size(elem,2)
                    f = @(xi,eta) function_cor(xi,eta).*phi{alpha}(xi,eta).*psi{beta}(xi,eta);
                    r = gauss_int2d_rec(f, V(1,:), V(2,:), V(3,:), V(4,:), N_gauss_int2d);
                    i = elem(n,beta);
                    j = elem(n,alpha);
                    matrix_A(i,j) = matrix_A(i,j) + r;
                end
            end
    end

end

end